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We have considered a one-dimensional coagulation-decoagulation system of classical particles on 
a finite lattice with reflecting boundaries. It is known that the system undergoes a phase transition 
from a high-density to a low-density phase. Using a matrix product approach we have obtained an 
exact expression for the average entropy production rate of the system in the thermodynamic limit. 
We have also performed a large deviation analysis for fluctuations of entropy production rate and 
particle current. It turns out that the characteristics of the kink in the large deviation function can 
be used to spot the phase transition point. We have found that for very weak driving field (when 
the system approaches to its equilibrium) and also for very strong driving field (when the system is 
in the low-density phase) the large deviation function for fluctuations of entropy production rate is 
almost parabolic while in the high-density phase it prominently deviates from Gaussian behavior. 
The validity of the Gallavotti-Cohen fluctuation relation for the large deviation function for particle 
current is also verified. 

PACS numbers: 05.40.-a,05.70.Ln,05.20.-y 

Keywords: entropy production, current fluctuations, large deviations in non-equilibrium systems, stochastic 
particle dynamics (theory) 



I. INTRODUCTION 

Physical systems in nature are either in or out of equi- 
librium. One-dimensional out of equilibrium systems 
usually exhibit unique critical behaviors such as phase 
transitions and shock formation [lj]. These properties 
make them very interesting to study from both math- 
ematical and physical point of views. Needless to say 
that many of these properties have not been fully under- 
stood yet. 

An out of equilibrium system is usually exposed to a flux 
of matter or energy. This results in a non-zero probabil- 
ity current between different microstates of the system 
and also entropy production. It is known that for a sys- 
tem in equilibrium, where detailed balance is hold, the 
entropy production is deterministically zero. This is the 
reason why the entropy production can be an indicator 
for being out of equilibrium. 

In recent years much attention has been paid to the stud 
of fluctuation theorems in out of equilibrium systems 
5]. It has been shown that for most of driven Markov 
processes the fluctuation theorem holds. Symmetry prop- 
erties of the large deviation functions of generalized cur- 
rent fluctuations have also been studied. These large 
deviation functions are also shown to obey a Gallavotti- 
Cohen (GC) type symmetry in systems with a finite state 
space. In contrast, for those systems with unbounded 
state space, such as the one-dimensional partially asym- 
metric zero-range process with open boundaries, the dis- 
tribution of large current fluctuations does not satisfy 
the GC symmetry ■ Other GC type symmetries have 
been found in a restricted class of Markov jump processes 
where the microscopic transitions have a particular struc- 
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ture and satisfy certain constraints [9(. 
In [l(| a large deviation analysis for fluctuations of par- 
tial and total particle currents in a zero-range process on 
a simple diamond lattice with open boundary conditions 
has been done. The validity of the GC fluctuation rela- 
tion for these particle currents is investigated and it has 
been found that the fluctuation relation is not satisfied 
for partial particle currents between sites even if it is sat- 
isfied for the total particle current flowing between the 
boundaries. 

So far, most of studies on the validity or breakdown of 
the GC symmetry and the fluctuation theorem have been 
mostly concentrated on boundary-driven systems. In this 
paper we consider an exactly solvable one-dimensional 
coagulation-decoagulation system. It is known that the 
system undergoes a phase transition in the steady state 
from a high-density to a low-density phase. We aim to 
study the large deviation functions for the fluctuations 
of entropy production rate and particle current and in- 
vestigate their symmetries in long-time limit. 
It turns out that the average entropy production rate 
of this system can be calculated exactly using a matrix 
product method. The average entropy production rate 
near the transition point changes abruptly, hence it can 
be used as a signal for the presence of a phase transition 
in the system. Studying the large deviation function for 
the entropy production rate reveals that the widely ob- 
served kink in the large deviation function disappears in 
the limit of large driving fields. The kink is also negligible 
as the driving field vanishes. However in the intermediate 
driving field regime the kink is prominent. This indicates 
that the characteristics of the kink in the large deviation 
function for the entropy production rate can be used to 
investigate the phase transition in the system. 
Our investigations also show that the large deviation 
functions for the local particle currents, defined as the 
net particle current through two consecutive lattice sites, 



2 



p= 1 -y 1 



I 1 1 1 

1 L 

FIG. 1: Simple sketch of a shock with two shock fronts on a 
lattice with L lattice sites. The density of particles in each 
region is denoted by p. 

do not satisfy the GC fluctuation relation. In contrast, 
the large deviation function for what we call the global 
particle current, defined as the sum of the local particle 
currents, satisfies the GC fluctuation relation. 
This paper is organized as follows: in the second sec- 
tion we will define the process and summarize the known 
results about its steady state. The second section is ded- 
icated to the mathematical basis of the fluctuation the- 
orem. Exact expression for average entropy production 
rate and its asymptotic behavior in the thermodynamic 
limit will be given in the fourth section. We will also 
perform a large deviation analysis for fluctuations of the 
entropy production in the fifth section. In the sixth sec- 
tion we will define the particle currents and check the 
validity of the GC fluctuation relation for these quanti- 
ties. We will finally summarize the results. 

II. DEFINITIONS AND KNOWN RESULTS 

Let us consider a system of interacting particles on a 
one-dimensional lattice of length L with reflecting bound- 
aries. The bulk of the system is assumed to evolve in time 
according to the following rules 
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in which A and stand for the presence of a particle and a 
hole in each lattice site respectively. Note that there is no 
injection or extraction at the boundaries. This particlc- 
nonconserving stochastic process has already been stud- 
ied extensively in related literature [TlT[l5| . The system 
has two different steady states: An empty lattice is a triv- 
ial steady state. It has been shown that the nontrivial 
steady state of this coagulation-decoagulation system is a 
matrix-product state (for a review of the matrix product 
approach see [16| ) provided that some constraints are ful- 
filled. More detailed investigations have shown that the 
steady state of the system can also be written in terms of 
a linear superposition of shock structures with two shock 
fronts (See FIG. Q] ). The shock fronts have simple ran- 
dom walk dynamics provided that some constraints on 



the microscopic reaction rates in ([1} are met. These re- 
sults confirm that the system has two different phases in 
the nontrivial steady state: a high-density (HD) phase 
and a low-density (LD) phase. In the HD phase the left 
(right) shock front moves preferentially to the left (right) 
while in the LD phase both shock fronts move preferen- 
tially to the left. 

For the sake of simplicity we adopt the following choices 
for the microscopic reaction rates in ([1} [If], [l2| 
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In terms of these newly defined microscopic reaction rates 
the system is in the HD phase for q < 7 while it is in the 
LD phase for q > 7. In the HD phase the bulk density of 
particles is equal to p = 1 — 7~ 2 while it is equal to zero 
in the LD phase. On the coexistence line where (7 = 7 the 
bulk density of particles on the lattice changes linearly. 
In the following section we will briefly review the math- 
ematical basis of fluctuation theorem. 



III. MATHEMATICAL PRELIMINARIES 

Let us consider a continuous-time Markov process with 
a configuration space which is denoted by 57. We assume 
that a spontaneous transition from configuration c to con- 
figuration c' takes place with a rate ui c _> c < where both c 
and c' belong to il. The time evolution of the probability 
distribution P(c,t), for the system being in c at time t, 
is given by a master equation which can be written as 

j t P(c,t) = -J2Hce'P(c',t), (3) 

where H is the Markov generator with elements 

{-Wc'^c if c ^ d 

(4) 

Using quantum Hamiltonian formalism [l[ the master 
equation ([3]) can be rewritten as 

l\P(t)) = -H\P(t)), (5) 

which is similar to the Schrodinger equation in imaginary 
time. 

The generating function for any current >J{t), which is 
a functional of the stochastic trajectory in configuration 
space f2, can be written as 

( e -^(*)) = ( S | e -^)*|P(0)), (6) 

in which (s| is a row vector with components (1, 1, 1, . . .) 
and that \P(0)) is the probability distribution vector at 
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t = 0. If the current J{t) changes its value by 9 C ^ C > 
whenever a jump from c —¥ d occurs, then the matrix 
elements of the modified generator H(n) in ([6]) can be 
written as B 



> c exp(-^f 



(7) 



In long-time limit the generating function 
ten as 
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in which e(/i) is given by the lowest eigenvalue of the 
modified generator ([7]). Finally the Legendre transfor- 
mation of e(/i), according to the Gratner-Ellis theorem, 
gives the large deviation function [j| 



e(x) = max„(e(/i) — x/x) . 



(9) 



According to the large deviation principle the probabil- 
ity distribution function, defined as the probability to 
observe a time-averaged value x = J(t)/t of the current 
J{t) over time interval [0,t], can now be written as 
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In limit of large time the probability distribution function 
V(x,t) satisfies 
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where E is a field conjugated to the flux J . The rela- 
tion (fTTj) can also be written as 



e(— x) — e(x) = Ex 



(12) 



which is also known as the CG fluctuation relation. De- 
pending on the definition of the time-integrated current 
J{t), the parameter x can be regarded as entropy pro- 
duction rate or particle current. 

Whenever the system jumps from c to d the entropy in 
the environment changes by In Wc ^ c ' [l7j . In order to 
calculate the total entropy changes in a trajectory one 
takes 8 c ^c' — In Wc ^ c ' in ([7]). The lowest eigenvalue of 
the modified generator ([7]) satisfies what we know as the 
GC symmetry 



e s (/i) = e s (l - fi) . 
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The large deviation function also satisfies the GC fluctu- 
ation relation 



e(-cr) - e(o-) = (S env )a, 



in which we have defined a := S/(S env )t. The average 
entropy production rate in the environment can be cal- 
culated using 
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In order to analyze fluctuations of particle current J in 
the steady state we need to count the net particle jumps 
during the time interval t. Whenever the system changes 
its configuration from c to d a particle might contribute 
to particle current. To count total particle jumps in a 
trajectory one can take 9 c ^,d — ±1 in © for those con- 
figuration changes which contribute to particle current on 
the lattice in two different directions [3J . In this case the 
lowest eigenvalue of the modified generator ([7j satisfies 
the GC symmetry 



ejG«) = ej{E - n) , 



(16) 



where E is the conjugate field. On the other hand, the 
large deviation function e(J) satisfies the GC fluctuation 
relation 



e(-J) - e(J) = EJ. 



(17) 



IV. AVERAGE ENTROPY PRODUCTION RATE 

If P*(c) is the probability of being in configuration c 
in the steady state, then besides the formula (fT5j) the 
average entropy production in the environment is given 
by 



c.c' 
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(18) 



provided that all transitions are reversible. Recent inves- 
tigations show that the behavior of the average entropy 
production in systems with out of equilibrium phase tran- 
sitions changes at the critical point [l8l - [23l |: therefore, it 
plays an important role in classifying different nonequi- 
librium phase transitions. 

In what follows we first show that (Tl8|) can be calculated 
exactly for the system defined by (JTJ and ([2]) . In order to 
calculate the steady state probability distribution P* (c) , 
one can use a matrix product approach. According to 
this approach P*{c) is written as a product of noncom- 
muting operators which satisfy an algebra [l6j . For the 
system defined by (TTJ) and © it has been shown that 
there exists a four-dimensional matrix representation for 
the quadratic algebra of the system [l2|, [l5[ . Using the 
four-dimensional matrix representation and the results 
obtained in [l^, [H| we find that in large-!/ limit the av- 
erage entropy production rate in the environment is given 
by 



( 2(<j 2 -l)(l-37 2 +2 7 4 )lng 



(14) (Senv) = < 



L 



for q < 7, 



for q > 7. 
(19) 

As can be seen (S env } changes discontinuously at the 
transition point. In other words, the average entropy pro- 
duction rate per lattice site in the steady state (S env )/L 
is a constant in the HD phase while it is zero in the LD 
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FIG. 2: Numerical results for a system of length L = 4: (a) The lowest eigenvalue of the modified generator for the entropy 
production for 7 = 2 and different values of q. (b) The average entropy production rate and its derivative for q = 2 as a 
function of 7. The inset shows its first derivative respect to 7. (c) The first derivative of the large deviation function for the 
entropy production for q = 2 and two different values of 7. In (d), (e) and (f) the large deviation function and its first and 
second derivatives are plotted for different values of q and 7 = 2. 



phase. This can be easily explained as follows: in the HD 
phase the lattice is nearly full of particles. This can be 
realized by the fact that in this phase the steady state is a 
linear superposition of product shock measures with two 
shock fronts in which left shock front moves preferentially 
to the left while the right shock front moves preferentially 
to the right (see FIG.Q]). The more particles contribute 
into the reactions, the more entropy is produced in the 
environment. Since nearly all lattice sites contribute in 
entropy production, the average entropy production (fT8|) 
in the HD phase is proportional to the system size. In 
contrast, in the LD phase both shock fronts move pref- 
erentially to the left (see FIG.[lJ; therefore, the lattice 
is almost empty. In this phase much less particles con- 
tribute in entropy production. That is why the average 
entropy production per lattice site is zero in large- L limit. 
We should also note that the average entropy production 
in the steady state (fl8|l is equal to zero for q = 1. This 
is the value of q for which the system is in equilibrium. 
We will discuss this later in forthcoming chapters. 



V. ENTROPY FLUCTUATIONS 

Finding an exact expression for the lowest eigenvalue of 
the modified generator for the entropy production e s (/i) 
of a system of length L is a formidable task; however, 
this can be done numerically for small lattices. Numeri- 
cally exact results obtained for a system of length L = 4 
are given in FIG. f2] We have plotted e s (/x) as a function 



of /1 for 7 = 2 in FIG. HJa). As q — > 1 the system ap- 
proaches to its equilibrium state, hence e s (/i) is almost 
parabolic [U, [2{| . It can be seen that e s (/i) behaves al- 
most the same way for large values of q. 
In FIG.[2jb) we have plotted the average entropy produc- 
tion rate obtained using (|15[) as a function of 7 for q = 2. 
The inset in this figure shows a peak in the derivative 
of the average entropy production rate which becomes 
more prominent by increasing system size L as (|19p has 
already predicted. This indicates that the behavior of 
(Senv) can be used to spot the transition point. 
In FIG. [He) the first derivative of large deviation func- 
tion for the entropy production rate respect to a for q = 2 
and two values of 7, one above and one below the tran- 
sition point, is plotted. Using (14]) one can easily find 
that HI 



e'{a) -e'(a) =2e'(a) + (S env ). 



(20) 



This shows that the sudden jump in the first derivative of 
the large deviation function is more pronounced for large 
values of the average entropy production rate i.e. in the 
HD phase. One should recall that the average entropy 
production rate as a function of 7 is of order L for q < 7 
while it is of order of unity for q > 7. 
In FIG. ©(d) we have plotted e(a) for 7 = 2 and dif- 
ferent values of q . It is known that the large devia- 
tion function for the entropy production rate exhibits a 
kink at a — which is a generic feature of the large 
deviation function and follows from the fluctuation the- 
orem [H [25| . As q — > 1 the kink at a — disappears 
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and the curve becomes a parabola which is, as we men- 
tioned, an indication for the system getting close to its 
equilibrium. It can be seen that e(cr) has almost the same 
behavior for q >> 1. As we will see in the next section 
the driving field which drives the system out of equilib- 
rium depend only on q. 

In FIG. d» and FIG. HJf) we have plotted the first and 
second derivatives of the large deviation function for the 
entropy production rate for different values of q. The ex- 
istence of a kink is best illustrated by these derivatives at 
a = 0. It can be seen that the jump in the first derivative 
of e(a) at a = disappears and its second derivative be- 
comes a constant as q — > I i.e. the system approaches an 
equilibrium steady state. The inset in FIG. [2jf) shows 
that as long as q =/= 1 the system is out of equilibrium. 
While being in the HD phase 1 < q < 7 the jump in the 
first derivative of the large deviation function increases 
as q increases. Note that (S env ) in the HD phase is an 
increasing function of q and of order of the system length 
L (see (HD)) an( i that the jump in the first derivative of 
the large deviation function is governed by {S env } as can 
be seen in (f^Uj) . 

In the LD phase q > 7 the average entropy production 
rate is negligible (of order of unity as can be seen in (JT9J)) 
and one expects that the jump in the first derivative of 
the large deviation function decreases in comparison to 
its value in the HD phase. This can also be seen in 
FIG. [2e) and FIG.^f) for q » 1. We expect that this 
will be more prominent as the system size is increased. 
We conclude this section by noting that the characteris- 
tics of the kink (the jump in first derivative of the large 
deviation function) can be used as a criteria for spotting 
the phase transition point. 

In the next section we will perform a large deviation anal- 
ysis for the particle current in the system. 



VI. PARTICLE CURRENT FLUCTUATIONS 

For the system defined by (fTJ) and @ the particle den- 
sity is not conserved. The time evolution of the av- 
erage local particle density (pk)(t) at a lattice site k 
(k = 1, ■ • • , L) is given by 



dt 



(p k )(t) = (J k ^)(t)-(J k )(t) + S k (t) (21) 



where {J k )(t) is called the average local particle current 
from lattice site k to k + 1 and S k is a source term. In 
the steady state the time dependency of the quantities 
will be dropped; therefore, the l.h.s of the equation (|2"Tj) 
becomes equal to zero and we find (see FIG. [3]) 



FIG. 3: Local particle currents and external sources in a sys- 
tem of length L. 



in which the average local particle current is defined as 

(Jk) = (q{p k pk+i) - q~ 1 (pkPk+i) 

+Ag((l - Pk)Pk+i) + - Pk)Pk+i) 

-Aq-^p^l - p k+1 )) - q-^Pkil - Pk+i))) 

= -q- 1 (l + A){p k )+q(l + A)(p k+1 ) 
-A(q - q~ 1 ){p k p k+1 ). 

It is easy to verify that Yl k =i $k = 0. We will also define 
an average global particle current as 

L-l 



(23) 



k=l 



Using the matrix product approach one can easily cal- 
culate the average local particle currents in the steady 
state. It turns out that the exact expression for the av- 
erage local particle current is given by 



(Jk) = 



'l - 7 2 ) (1 - q 2 ) 7 -2fc+2i-4 g -4fc+2i-3 



S k = (Jk) - (Jk-i) for k = 1, 



,L. 



(22) 



7 2L + (jq) 2L (q 2L - 7 2L ) - q 2L 

x ( 7 2fe 9 2L (7 2 (7 2 + 9 2 -l)(79) 2 ' £ 
+ 7 2 (l-7 2 )(l + (? 4 ) _ g 2 ) 

+ q 2k 1 2L {q 2 {l-2 1 2 )( 1 q) 2k 

+ 7V(7V-9 2 + l))) (24) 

for k = 1, • ■ • , L - 1. In FIG. H we have plotted flM]) as 
a function of 7 for q — 2 and L = 100. It can be seen 
that for 7 < 2 i.e. in the LD phase, the average local 
particle currents are zero except in the vicinity of the left 
boundary. However, for 7 > 2 i.e. in the HD phase, the 
average local particle currents are nonzero throughout 
the lattice. 

In order to calculate the average local particle currents 
one can equivalently construct a modified generator for a 
jump process which counts the number of local or global 
particle jumps on a trajectory over a time interval [0, t]. 
This can be done using ([7]) and an appropriate choice 
for # c _>. c < as we explained in mathematical preliminaries 
section. The first derivative of the lowest eigenvalue of 
this operator respect to p at p, = gives the average 
particle current. 

We have found that only the global particle current J 
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FIG. 4: The average local particle currents ((24} as a function 
of 7 for q = 2 and L = 100. 

satisfies GC fluctuation relation (|17p and that its average 
is proportional to the average entropy production rate 

(S mv ) = E(J) (25) 

where E is equal to lng 2 . The thermodynamic force E 
which is conjugated to the flux of particle (J) physically 
corresponds to the effective driving field pushing parti- 
cles. Our numerically verified results show that the low- 
est eigenvalue of the modified generator for the global 
particle current and that of the entropy production sat- 
isfy 

e,(ji) = ej{E^). (26) 

The fact that the local particle currents do not fulfill 
GC fluctuation relation has already been observed in 
a zero-ran ge p rocess on a diamond lattice with open 
boundaries [T(j- In this paper we study an exactly solv- 
able system with a particle-nonconserving dynamics. 
Using the fact that the system is at equilibrium when 
q = 1 we can calculate the conjugate field E as follows. 
We imagine that there exists an equilibrium system with 
rates w e ^ c , which obey detailed balance. Now by apply- 
ing an external field E, we recover our system defined 
in (JTJ and @ whose microscopic transition rates satisfy 
the following relation 

«Wc = w^eT^ (27) 

in which 6 c ^ c i — 1 {9 r ^. c > = —1) if the transition from c 
to d is associated with a particle jump to the left (right) 
on the lattice. The relation (|27p gives E = In q 2 provided 
that the microscopic rates of the system in equilibrium 
are given by Q with q = 1. 

The matrix product approach predicts that the equilib- 
rium probability distribution function P* q (c), i.e. the 
probability distribution function P*(c) at q = 1, has the 
following properties 

p* g (---io---) = i 5 *0--oi---), 

P e * g (---ll---) = (7 2 -l)^(---01---)- (28) 



ft is not hard to verify that the equilibrium probability 
distribution function satisfies the local detailed balance 
condition given by Q 

uwe-f e — 'P* 9 (c) = w ct ^ c e-^'^p* q (d). (29) 

It is known that the local detailed balance condition leads 
to the GC symmetry of the global particle current [2j. 
Let us now have a look at the behavior of (J) as a func- 
tion of q. As we mentioned above, the driving force E 
is zero at q = 1 and therefore the system is at equilib- 
rium. At this point (J) is zero. As q is increased the 
driving force increases and hence the system is driven 
out of equilibrium. The average global particle current 
is an increasing function of q for 1 < q < 7 i.e. in the 
HD phase. The phase transition occurs at q = 7. Above 
the transition point i.e. in the LD phase, (J) becomes 
negligible. 

In FIG. [5] we have plotted the large deviation function 
for the global particle current and its first derivative for 
q = 2 and three values of 7 in a system of length L = 4. 
These point are chosen to be above, below and at the 
transition point. As can be seen, a sudden jump exists 
in the first derivative of the large deviation function at 
J = 0. It can be seen that in the LD phase the minimum 
of the large deviation function occurs at a point (which 
gives the average global particle current) close to zero 
while it is nonzero in the HD phase. This is in accor- 
dance with the average entropy production rate behavior 
as we explained before. 

In order to check the validity of the GC fluctuation re- 
lation for the global and local particle currents we have 
plotted (|17p both for J and J^s in a system of length 
L — 4 in FIG. \5\ It can be seen that only the global 
particle current satisfies the GC fluctuation relation in 
accordance with (fl~7|) . The slop of this line is equal to E. 
For the local particle currents a linear behavior can be 
seen only in the vicinity of the origin. 

VII. CONCLUSION 

In this paper we have considered a one-dimensional 
classical system with reflecting boundaries and a particle- 
nonconserving dynamics. It is known that by varying the 
microscopic reaction rates, the system undergoes a phase 
transition from a LD phase to a HD phase. 
Using a matrix product approach we have obtained exact 
expression for the average entropy production rate in the 
environment in the long time limit. It turns out that the 
average entropy production rate changes discontinuously 
at the phase transition point while it is zero at q = 1 
where the system is in equilibrium. 

We have studied the entropy fluctuations in the system 
for L = 4. We have found that the large deviation func- 
tion for the entropy production rate becomes a parabola 
as the system approaches to its equilibrium. At zero 
entropy production rate a kink is observed in the large 
deviation function for the entropy production rate. The 
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FIG. 5: (Left to right) The large deviation function for the global particle current and its first derivative in a system of length 
L = 4 for q = 2 and three values of 7. These points are located below, above and at the transition point. The last figure shows 
that the large deviation function for the global particle current satisfies the GC fluctuation relation theorem (|17[1 . 



kink disappears at both very large and very small driving 
fields. We expect that in large- L limit the kink (also the 
discontinuity in the first derivative of the large deviation 
function for the entropy production rate) become more 
noticeable in the HD phase. 

We have also investigated the validity of the GC fluc- 
tuation relation for the particle current in this system. 
We have considered two types of particle currents: local 
particle currents which are defined as particle currents 
between consecutive lattice sites and global particle cur- 
rent as a sum of these local particle currents. The average 
particle currents are calculated exactly. Our numerical 
investigations reveal that only the global particle current 
fulfills the GC fluctuation relation. Moreover, our analyt- 



ical and numerical investigations show that the average 
global particle current (J) is proportional to the average 
entropy production rate (S env ). 

Previous investigations have shown that the first deriva- 
tive of the average entropy production rate displays a 
peak, a discontinuity or a divergence at criticality. In 
this paper we have shown that the stationary average 
entropy production rate in our system, defined by ([T]) 
and ((2J, changes discontinuously at the critical point. It 
seems that different nonequilibrium phase transitions can 
be classified using criticality of average entropy produc- 
tion rate at a transition point. 
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